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The finest state space resolution that can be achieved in a physical dynamical system is limited by 
the presence of noise. In the weak-noise approximation the stochastic neighborhoods of deterministic 
periodic orbits can be computed from distributions stationary under the action of a local Fokker- 
Planck operator and its adjoint. We derive explicit formulae for widths of these distributions in the 
case of chaotic dynamics, when the periodic orbits are hyperbolic. The resulting neighborhoods form 
a basis for functions on the attractor. The global stationary distribution, needed for calculation of 
long-time expectation values of observables, can be expressed in this basis. 
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I. INTRODUCTION 

This paper investigates the interplay of deterministic 
chaotic dynamics and weak stochastic noise, and pro¬ 
poses a new definition of the neighborhood of a noisy 
hyperbolic state space point. Such neighborhoods are 
conjectured to partition the state space in optimal way 
and provide a basis function set for the evaluation of the 
stationary distribution. 


A. Width of a noisy trajectory 

The basic idea of a stochastic ‘neighborhood’ is that 
the balance between the noise broadening of a trajec¬ 
tory and the deterministic contraction leads to a proba¬ 
bility distribution of finite width, as opposed to one that 
spreads with time (diffusion only). For an orbit that con¬ 
verges to a linearly stable, attractive equilibrium, this 
neighborhood was computed in 1810 by Laplace [ , 2] 
and is today known as a solution to the Lyapunov equa¬ 
tion [3], or the Ornstein-Uhlenbeck process [4]: for a 1- 
dimensional flow, the deterministic equilibrium point is 
smeared into a Gaussian probability density centered on 
it, whose covariance Q = —D/A is a balance of the expan¬ 
sion rate D (diffusion constant) against the contraction 
rate A < 0. Fokker-Planck equation [5] generalizations to 
higher-dimensional stable equilibria and limit cycles (sta¬ 
ble periodic orbits) are immediate, provided proper care 
is taken of the diffusion along the periodic orbit [6, 7]. 

What if a periodic orbit is unstable? Both the diffusion 
rate and the linearized stability rate A > 0 now expand 
forward in time, and cannot balance each other. This 
problem was solved in refs. [8, 9] for repelling periodic 
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orbits with no contracting directions , by balancing the 
stochastic diffusion against the contraction by the adjoint 
Fokker-Planck operator. The resulting covariance matrix 
defines the stochastic neighborhood for a repelling orbit, 
while the Ornstein-Uhlenbeck covariance defines it for a 
stable orbit. However, neither these stable nor repelling 
orbits play a role in chaotic dynamics. The long-time 
attractors of chaotic dynamics are organized by an in¬ 
finity of hyperbolic periodic orbits [ 0—12] , orbits which 
an ergodic trajectory visits by approaching them along 
their stable eigendirections, and leaves along their unsta¬ 
ble eigendirections. 

The central result of this paper is that techniques de- 
velopped for solving the Lyapunov equations [ 3-1 ] en¬ 
able us to define the neighborhood of a hyperbolic peri¬ 
odic point by splitting the covariance matrix Q into two 
(mutually non-orthogonal) covariance matrices, Q cc for 
contracting directions, and Q ee for the expanding direc¬ 
tions. 

There are two immediate applications of the notion 
of the neighborhood of a hyperbolic point: (a) ‘optimal 
partition’ of the attractor, and (b) construction of a basis 
set for the stationary distribution of a noisy chaotic flow. 


B. An optimal partition from periodic orbits 

While in the idealized deterministic dynamics the state 
space can be resolved arbitrarily finely, in physical sys¬ 
tems noise always limits the attainable state space reso¬ 
lution. 

This observation had motivated the many limiting res¬ 
olution estimates for state space granularity of chaotic 
systems with background noise. The idea of an opti¬ 
mal partition in this context was first introduced in 1983 
by Crutchfield and Packard [16] who formulated a state 
space resolution criterion in terms of a globally averaged 
“attainable information.” The approach was later gen¬ 
eralized and applied to time-series analysis, where the 
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underlying dynamics is unknown [17, 18]). A different 
strategy consists of computing a transfer matrix between 
intervals of a uniform grid, and estimating averages of ob¬ 
servables from its eigenvalues and eigenfunctions. First 
introduced by Ulam [ 9], this technique has been devel¬ 
oped over the years [20, 21]. All of these approaches (see 
ref. [9] for a review) are based on global averages, and as¬ 
sume that granularity is uniform across the state space. 
In contrast, the main, computationally precise lesson of 
our work is that even when the external noise is white, 
additive, and globally homogenous, the interplay of noise 
and nonlinear dynamics always results in a local stochas¬ 
tic neighborhood, whose covariance depends on both the 
past and the future noise integrated and non-linearly con¬ 
volved with deterministic evolution along the trajectory. 
The optimal resolution thus varies from neighborhood to 
neighborhood, and has to be computed locally. As was 
shown in ref. [ ] for a strictly expanding 1 d chaotic map 
and a given noise, the maximal set of non-overlapping 
neighborhoods of periodic orbits can be used to construct 
an L optimal partition ’ of the state space, and compute dy¬ 
namical averages from the associated approximate matrix 
Fokker-Planck operator. 


C. The stationary distribution 

In this paper we utilize our construction of optimal 
partitions to approximate the stationary probability dis¬ 
tribution function by a finite sum over Gaussians, one for 
each neighborhood. When the dynamics is chaotic, the 
most one can predict accurately for long times are the 
statistical properties of the system, given by the state- 
space averages of observables a(x), 

<a> = / dxp(x)a(x ), (1) 

where the stationary distribution (natural measure [22- 
24]) p{x) is the probability of finding the system in the 
state x on the attractor. For a deterministic system p(x) 
is a singular, nowhere differentiable distribution with 
support on a fractal set, and its numerical computation 
is usually not feasible. However, for any physical system 
the noise washes out fine details of the dynamics, and 
the stationary distribution is smooth. Here we propose 
a smooth function basis for p(x), based on the optimal 
partition of the state space. We develop our formal¬ 
ism for discrete-time dynamical systems and illustrate 
it by computing the neighborhoods and estimating p(x) 
for the Lozi map [25], a simple 2-dimensional discrete¬ 
time chaotic system. The idea is to first partition the 
attractor into an optimal partition set of neighborhoods, 
and then use the associated local Gaussian distributions 
as a finite set of basis functions for the global station¬ 
ary distribution. In the 2-dimensional Lozi example, our 
estimates for the stationary distribution are consistent 
with those obtained by the direct numerical estimation 


of the lattice-discretized probability densities computed 
from long stochastic (Langevin) trajectories. 


II. THE NEIGHBORHOOD OF A HYPERBOLIC 
POINT 

An autonomous discrete-time stochastic dynamical 
system (AJ,/, A) can be defined by specifying a state 
space M, a deterministic map / : A4 —> AJ, and an addi¬ 
tive noise covariance matrix (diffusion tensor) A = A(x). 
In one time step, an initial Dirac-delta density distribu¬ 
tion p a (x) located at x is smeared out into a Gaussian 
ellipsoid p a +i(y) centered at y = /(x), with covariance 
A(x). This defines the kernel of the Fokker-Planck evo¬ 
lution operator in d dimensions [5] 

C FP (y,x) dx = e - 2 a(*) ( y~f ( x )) 

[dx] = d d x/ det(27rA) 1 / 2 . (2) 

Consider a trajectory {x a } = (x a ,x a +i,x a +2, * * ■) gener¬ 
ated by the deterministic evolution rule x a+ i = f(x a ), 
and shift the coordinates in each x a neighborhood to 
x = x a + z a . In the vicinity of x a the dynamics can 
be linearized as z a +1 = M a z a , where M a = df(x a ) is 
the one time-step Jacobian matrix. 

Prepare the initial density of trajectories p a (z a ) in the 
x a neighborhood as a normalized Gaussian distribution 
p(z a ,Q a ), centered at z a = 0, with a strictly positive- 
definite covariance matrix Q a . The support of density 
p(zaiQa) can be visualized as an ellipsoid with axes 
oriented along the eigenvectors of Q a . The linearized 
Fokker-Planck operator 

C(z a+ 1 ,Z a )dz a = e -(^ +i -M aZa y^( Za+t -M aZ a) [dZa] 

maps this distribution one step forward in time into an¬ 
other Gaussian 

P(Za-\-l ? Qa-\-l) = J dz a C, (z a _|_i, Z(f) p(z a: Q a ) , (3) 

with the covariance matrix deformed by the dynamics 
and spread out by the noise, as given by the discrete 
Lyapunov equation [3, 26], 

Qa+1 = MaQaM a T + A a . (4) 

In other words, the two covariance matrices, (i) the de¬ 
terministically transported Q a —)> M a Q a M a T , and (ii) 
the noise diffusion tensor A a , add together in the usual 
manner, as squares of errors. 

Similarly, density evolution for dynamics with strictly 
expanding Jacobian matrices M a can be described by the 
action of the adjoint Fokker-Planck operator [8, 9], with 
kernel 

C\y, x) dy = e ~ 2 (y~f( x )) ar( y~f ( x )) _ 
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The adjoint Fokker-Planck operator expresses the current 
density p a as the convolution of its image p a +\ with the 
noisy dynamics 


Pa^aiQa) J dz a -\-\£^ (z a •> ^a+l )Pa+l (^a+1 5 Qa-\-l ) • 

Like in the forward evolution, we may substitute a Gaus¬ 
sian density into this equation to obtain the discrete ad¬ 
joint Lyapunov equation for the covariance matrices 

M a QaMa T = Qa +1 + A a • (5) 

We show in what follows that, if the Jacobian matrices 
M a have all eigenvalues strictly contracting (expanding), 
any initial Gaussian converges to an invariant density 
under the action of the (adjoint) Fokker-Planck operator. 
Consider first the case of a map f(x) with a stable fixed 
point at x a (at z = z a = 0). The covariance matrix 
transforms as 

Q = A + MAM T + M 2 A(M T ) 2 + ■■■ 

oo 

= Y, M n A(M T ) m 5 mn . (6) 

m,n =0 

By inserting the Fourier representation of Kronecker J mn 
into (6), we can recast this expression into the resolvent 
form 


Here the blocks A e and A c contain all expanding, con¬ 
tracting eigenvalues of the monodromy matrix, respec¬ 
tively. The covariance matrix Q = SQS T is not block- 
diagonalized by the above similarity transformation, but 
consider the four blocks 


A - f d0 1 A 1 

Qij -J 0 2tt 1 - e-*Ai ij 1 - e*Aj ’ 


( 10 ) 


where A = S~ 1 A(S~ 1 ) T , and i, j G {c, e}, where {c, e} 
denotes {contracting, expanding}. This expression may 
be evaluated as a contour integral around the unit circle 
in the complex plane [13, 15] 

® ij = 2^f v i - zaj ■ (11) 


The diagonal blocks Q cc , Q ee have either all expand¬ 
ing or all contracting eigenvalues, meaning at least one 
pole inside and one pole outside the unit circle, and the 
residue theorem yields a non-vanishing result for the in¬ 
tegral. Consider next the off-diagonal block Q ce withA^ 
contracting and A j expanding: in this case the poles all 
lie outside the unit circle, and the integral vanishes. The 
remaining off-diagonal block having A i expanding and A j 
contracting must also vanish when integrated, due to the 
symmetry of Q, which is therefore block-diagonal, 


n2n 7 n °° 

Q= — E (e~ ie M) n A(e ie M T y 

'' 0 m n—n 


m,n=0 

1 


t 2v d0 

J 0 2-rr 1 - e- ie M^l - e ie M T ' 


( 7 ) 


We do the same in the expanding case, by using the ad¬ 
joint evolution 

o° i -j 

MQM t = V —A— pp— 5 nm 
c M n (M T ) m 

m,n=0 v 7 


-l 

-L 


0 27r 
27T 


dd_ y, 

m,n =0 

de i 




o —i0 

\MT 

1 


o 2ir 1 — e lS /M 1 — e~ ie /M J 


, ( 8 ) 


Q = s ( Q o l„) sT - (12> 

These results are easily extended to a periodic orbit p 
of period n p , since any point x a of the orbit is a fixed 
point of the n p th iterate of the map. The forward and 
adjoint evolution equations (4) and (5) for the covariance 
matrix, as well as the resolvent (7) all still hold, with 
some changes in the notation: each periodic point x a has 
its own neighborhood, with its own covariance matrix 
Q a . The monodromy matrix M a of x a now evolves n p 
steps along the orbit 

Ma P — M a +n p -1 • ' • Af a+ 2Af a+ iM a , 

while the diffusion tensor A a now accounts for the total 
noise accumulated along the periodic orbit, 


which is then easily reduced to (7), so that the resolvent 
form is the same regardless of whether M is expanding or 
contracting. This result comes particularly handy when 
we deal with a hyperbolic fixed point, that is when M a 
has both expanding and contracting eigenvalues. The 
monodromy matrix is not symmetric, and it cannot be 
diagonalized by an orthogonal transformation, but its ex¬ 
panding and the contracting parts can be separated with 
a similarity transformation S that brings M to a block- 
diagonal form, 

A S S-'MS = ( A 0 > l ) . (9) 


A p,a = E K+i+1 1 A -+i K+i+1 1T - ( 13 ) 

i =0 

III. OPTIMAL PARTITION AND STATIONARY 
DISTRIBUTION 

At this point our strategy is to build a partition out 
of neighborhoods of the periodic points, each defined by 
means of the stationarity condition (7): solve for the ex¬ 
panding and contracting blocks of (12) separately, and 
draw a parallelogram on the supports of the resulting 
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FIG. 1: Building the partition for the Lozi attractor, for an 
isotropic constant diffusion tensor Sij A. In frames (a) and (b) 
A = 0.01. The deterministic attractor is the fractal structure 
in the background of each picture. Stochastic neighborhoods 
of a set of periodic points are indicated by their standard 
deviation parallelograms, (a) An initial partition: with only 
periodic points of periods < 5 much of the attractor remains 
to be covered, (b) The final, optimal partition covers the 
whole attractor, with no pair of neighborhoods overlapping by 
more than 50%. (c) N, the number of neighborhoods needed 
to achieve the optimal partition for a given noise strength A. 


Gaussians, with axes oriented along the eigenvectors of 
the covariance matrices Q ee and Q cc , their widths given 
by one standard deviation along each direction. We say 
that two neighborhoods overlap if they do so by at least 
50% of their areas (consistent with the la confidence in¬ 
terval chosen as overlapping threshold in ref. [ ]). For 
a typical chaotic map, periodic points are dense in the 
deterministic attractor [27], which we now aim to cover 
entirely with the minimum number of neighborhoods pos¬ 
sible. We do so via the following algorithm: i) Find 
periodic points of period n p = 1 , 2 ,..., and their corre¬ 
sponding neighborhoods, ii) If any neighborhood over¬ 
laps with the neighborhood of a shorter periodic point, 
then it is discarded and the neighborhood of lower period 
occupying the same area is instead kept in the partition. 
Hi) Among groups of neighborhoods of the same period, 
discard those that overlap, while keep the rest in the 
partition, iv) The algorithm stops when the attractor is 
fully covered and no further non-overlapping neighbor¬ 
hoods can be found. An example is shown in Fig. 1 for 
the two-dimensional Lozi attractor [25]. 

The main utility of a good partition is that it provides 
a basis for an accurate and efficient estimate of long-time 
averages of observables defined on the dynamical system, 


of the form ( 1 ). As explained in sect. IC, our goal here is 
to determine the stationary distribution p{pc). For that 
purpose, we use as basis Gaussian ellipsoids that satisfy 
the local stationarity condition (7) in each neighborhood 
of the optimal partition. A set of Gaussians centered 
at every point in the state space M forms an overcom¬ 
plete, non-orthogonal basis for functions in L 2 (J\A), as is 
well known from the study of coherent states of quan¬ 
tum harmonic oscillators [28]. Our (also overcomplete 
and non-orthogonal) set of Gaussians is centered only on 
periodic points, which are dense in the deterministic at¬ 
tractor, but not in the entire state space. Therefore our 
basis is designed to resolve the structure of any function 
with support on the hyperbolic ‘strange set’ (an attrac¬ 
tor or a repeller). The Gaussians are constructed so that 
their widths balance the noise spreading and the (time- 
forward or -backward) contraction of the deterministic 
dynamics. In the transverse directions, the basis gives 
the width of the global stationary distribution, locally 
everywhere determined by the balance between noise and 
dynamics. Along the attractor, the basis determines the 
minimum number of neighborhoods needed to fully re¬ 
solve the structure of the stationary distribution. 

There are numerical methods (such as refinements of 
Ulam’s method [29]) that identify the asymptotic attrac¬ 
tor by running long noisy trajectories, dropping the tran¬ 
sients, and covering the attractor so revealed by a finite 
number of boxes. These algorithms have no a priori in¬ 
formation about how the stationary distribution behaves 
transversely to the deterministic attractor, and they may 
easily overestimate the number of basis elements needed 
to resolve this structure. In contrast, in our approach 
the transverse structure is automatically accounted for 
by the local balance between the noise and the determin¬ 
istic contraction along the stable, transverse directions, 
given by covariance matrix block Q cc in ( 12 ). Further¬ 
more, estimating p{pc) by binning long noisy trajectory 
over a finite number of attractor-covering boxes is feasi¬ 
ble only in a low-dimensional state space, while ( 12 ) can 
be computed for state space of any dimension. 

In discrete time dynamics, the stationary distribution 
is the ground-state eigenfunction of the Fokker-Planck 
evolution operator ( 2 ) with escape rate 7 , 

C p{pc) = e~ 7 p(x) . (14) 

In order to estimate the stationary distribution, we write 
it as a sum over the neighborhoods of the periodic points: 

N 

P{ x ) = ^2 h a<t>a{x), (15) 

a= 1 

where </> a = e~ x * ® aXa are the Gaussian basis functions, 
with Q a given by ( 12 ), and the coefficients {h a } to be 
determined. The truncation of the expansion (15) to N 
basis functions follow from our optimal partition. We es- 
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FIG. 2: (Color online) The stationary distribution of the 
Lozi map with A = 0.01. (a) A direct numerical calcula¬ 

tion obtained by running noisy trajectories for a long time, 
(b) The stationary distribution calculated with the optimal 
partition method. 


timate the coefficients h a by minimizing the cost function 



2 

dx , 


(16) 


together with the normalization constraint for p(x). We 
can also estimate the escape rate of the system by mini¬ 
mizing the error with respect to e -7 . 

As an example, we apply the procedure to the Lozi 
map [2 ] 


^n+l — 1 ^|^n| T 6 Vn 

Vn+l = %n (17) 


with parameters a = 1.85,6 = 0.3 and isotropic, con¬ 
stant diffusion tensor Sij A with A ranging in the inter¬ 
val [0.003,0.1]. Figure 1 (c), which shows the number N 
of neighborhoods required by the optimal partition for a 
given A, illustrates the efficiency of our method: N goes 
from tens to few hundreds in the noise range considered. 
In order to test our algorithm, we also estimate p(x) and 
7 by a direct numerical simulation. The (x, y) state space 
is divided into uniform mesh 6.4 x 10 5 bins; we follow long 
stochastic trajectories and count how many times they 
visit each bin. The stationary distribution ps{x) is then 
the normalized frequency distribution of the whole grid. 
The deterministic Lozi map has a fixed point at the edge 
of attractor, whose stable manifold is the boundary of 
the deterministic basin of attraction. The noise makes it 
possible for a stochastic trajectory to cross this boundary 
and escape. We compute the escape rate as the ratio of 
the total number of points in the noisy trajectories to the 
number of escapes. Fig. 2 shows an example of the sta¬ 
tionary distribution estimated with both methods, while 
in Fig. 3 we compare quantitatively the two procedures. 
In particular, we estimate the relative error between the 
stationary distribution p computed with the optimal par¬ 
tition and pB computed on the uniform grid, by using a 
normalized L 2 distance, as 


d(p, pb) 


f (p(x) ^PBpJpfdx 

f (Pb(x )) 2 dx 


(18) 



FIG. 3: (a) The escape rate from the attractor as the function 
of the noise strength A. Squares: uniform grid discretization 
method. Triangles: optimal partition, (b) The normalized L 2 
distance d(p, ps) between the corresponding stationary distri¬ 
butions. 


The two distributions are within 5% of each other, 
whereas the escape rates differ by at most 10 %, over a 
range of A that spans two orders of magnitude. 


IV. DISCUSSION 


In conclusion, we have generalized the optimal parti¬ 
tion hypothesis first formulated in [8] to hyperbolic maps 
in arbitrary dimension, and tested the method on a two- 
dimensional system with weak white isotropic noise. As 
noise induces a finite resolution of the state space of any 
physical system, finite numbers of neighborhoods suffice 
to partition the state space explored by chaotic dynam¬ 
ics, and to estimate long-time averages of observables. 
Here we have used the deterministic unstable periodic 
orbits as the skeleton on which to build an optimal par¬ 
tition for the noisy state space. First we determine a 
local stationary distribution in the neighborhood of each 
periodic point by balancing the noise against the deter¬ 
ministic expansion or contraction. From the separation 
of expanding and contracting blocks in the covariance 
matrix that characterizes the Gaussian approximation to 
the local stationary distribution, we carve out a precise 
definition of neighborhood, the constituent of our parti¬ 
tion, which is then used to approximate the global sta¬ 
tionary distribution, estimate the escape rate (for open 
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systems that allow escape), and any long-time averaged 
observable. Numerical tests confirm that the accuracy of 
our method is comparable to that of a uniform grid dis¬ 
cretization, but the number neighborhoods required for 
our optimal partition 10 to 100) is three-four orders 
of magnitude smaller than the number of bins used in the 
uniform grid discretization method 10 5 ). 

The problems dynamical chaos (or ‘turbulence’) the¬ 
ory faces nowadays are not two- but high-, even infinite¬ 
dimensional. Today it is possible to compute numerically 
exact periodic orbits (‘recurrent flows’ [30]) in a variety 
of physically realistic turbulent fluid flows [31, 32], but 
these calculations are at the limit of what current codes 
can do, and we hope that the methods presented here 
can provide sharp criteria for when a sufficient number of 
such solutions has been computed. Furthermore, unlike 
the uniform grid discretization, our partitions are smart, 
since they rely on the periodic orbits of the determinis¬ 
tic system as skeleton of the dynamics, as well as effi¬ 
cient, due to the finite (and optimal!) numbers of neigh¬ 
borhoods and corresponding basis functions. This, we 
believe, should make our algorithm less costly to imple¬ 
ment than direct numerical simulations in higher dimen¬ 
sions, where discretizations would be impractical. With 
some modifications and application of Poincare sections, 


the formalism can be applied to continuous time flows as 
well [9, 33]. Outstanding challenges include dealing with 
the lack of hyperbolicity in higher dimensions (marginal 
directions were treated in ref. [9] for Id maps), as well as 
extending the definition of neighborhood to other time- 
invariant sets, such as relative periodic orbits and par¬ 
tially hyperbolic invariant manifolds. Further technical 
issues, such as improving the efficiency of the minimiza¬ 
tion algorithm by modifying the basis of functions used in 
the computation of the stationary distribution, are also 
part of our agenda. 
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